clear all;clc;
tic

er0=8.854187817e-12;
mu0=4*pi*1e-7;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NodePt=dlmread('node.txt');
PlaneNd=dlmread('plane.txt');
SurfaceNd=dlmread('surface.txt');

GlobKe=zeros(length(NodePt));
GlobBe=zeros(length(NodePt),1);

J=1e6;
alfa=1;

NodeNum=size(NodePt);
PlaneSize=size(PlaneNd);

figure
subplot(1,3,1)
hold on
axis equal
for n=1:PlaneSize(1);
    for np=1:3
        nodex(np)=NodePt(PlaneNd(n,np+2),2);
        nodey(np)=NodePt(PlaneNd(n,np+2),3);
    end
   
    be(1)=nodey(2)-nodey(3);
    be(2)=nodey(3)-nodey(1);
    be(3)=nodey(1)-nodey(2);
    
    ce(1)=-nodex(2)+nodex(3);
    ce(2)=-nodex(3)+nodex(1);
    ce(3)=-nodex(1)+nodex(2);
      
    plot([nodex(1) nodex(2) nodex(3) nodex(1)],[nodey(1) nodey(2) nodey(3) nodey(1)])
    
    deltae=(be(1)*ce(2)-be(2)*ce(1))/2;
    
    Keij=zeros(3,3);
    for ni=1:3;
        for nj=1:3
            Keij(ni,nj)=alfa*(be(ni)*be(nj)+ce(ni)*ce(nj))/(4*deltae);
            GlobKe(PlaneNd(n,ni+2),PlaneNd(n,nj+2))=GlobKe(PlaneNd(n,ni+2),PlaneNd(n,nj+2))+Keij(ni,nj);
        end
    end
    
    bei=zeros(1,3);
    if(PlaneNd(n,2)==1)
        for ni=1:3
            bei(ni)=mu0*deltae*J/3;
            GlobBe(PlaneNd(n,ni+2))=GlobBe(PlaneNd(n,ni+2))+bei(ni);
        end
    end
end

SurfaceSize=size(SurfaceNd);
NdFlag=zeros(1,NodeNum(1));

for n=1:SurfaceSize(1) 
    for np=1:2
        nodex(np)=NodePt(SurfaceNd(n,np+1),2);
        nodey(np)=NodePt(SurfaceNd(n,np+1),3);
    end
    plot([nodex(1) nodex(2) nodex(1)],[nodey(1) nodey(2) nodey(1)],'r','linewidth',2)
     
    if(NdFlag(SurfaceNd(n,2))==0)
    GlobKe(SurfaceNd(n,2),SurfaceNd(n,2))=1e15;
    GlobBe(SurfaceNd(n,2))=0;
    NdFlag(SurfaceNd(n,2))=1;
    end
    
    if(NdFlag(SurfaceNd(n,3))==0)
    GlobKe(SurfaceNd(n,3),SurfaceNd(n,3))=1e15;
    GlobBe(SurfaceNd(n,3))=0;
    NdFlag(SurfaceNd(n,3))=1;
    end
end
NdFlag=zeros(NodeNum(1),1);

phi=GlobKe\GlobBe;
X=NodePt(:,2);
Y=NodePt(:,3);
subplot(1,3,2)
tri0=delaunay(X,Y);
trisurf(tri0,X,Y,phi)
axis equal

Bx=zeros(length(NodePt),1);
By=zeros(length(NodePt),1);
Ex=zeros(length(NodePt),1);
Ey=zeros(length(NodePt),1);

for n=1:length(PlaneNd);
    for np=1:3
        nodex(np)=NodePt(PlaneNd(n,np+2),2);
        nodey(np)=NodePt(PlaneNd(n,np+2),3);
    end
   
    be(1)=nodey(2)-nodey(3);
    be(2)=nodey(3)-nodey(1);
    be(3)=nodey(1)-nodey(2);
    
    ce(1)=-nodex(2)+nodex(3);
    ce(2)=-nodex(3)+nodex(1);
    ce(3)=-nodex(1)+nodex(2);
    
    deltae=(be(1)*ce(2)-be(2)*ce(1))/2;
    
    for ni=1:3
        for np=1:3
            Bx(PlaneNd(n,2+ni))=Bx(PlaneNd(n,2+ni))+ce(np)*phi(PlaneNd(n,np+2))/2/deltae;
            By(PlaneNd(n,2+ni))=By(PlaneNd(n,2+ni))-be(np)*phi(PlaneNd(n,np+2))/2/deltae;

            Ex(PlaneNd(n,2+ni))=Ex(PlaneNd(n,2+ni))+be(np)*phi(PlaneNd(n,np+2))/2/deltae;
            Ey(PlaneNd(n,2+ni))=Ey(PlaneNd(n,2+ni))+ce(np)*phi(PlaneNd(n,np+2))/2/deltae; 
        end
        NdFlag(PlaneNd(n,2+ni))=NdFlag(PlaneNd(n,2+ni))+1;
    end
end

Bx=Bx./NdFlag;
By=By./NdFlag;
Ex=Ex./NdFlag;
Ey=Ey./NdFlag;

subplot(1,3,3)
hold on
quiver(X,Y,Bx,By)
axis equal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
toc